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ABSTRACT 



Context. The evolution of dust particles in protoplanetary disks determines many observable and structural properties 
of the disk such as the spectral energy distribution (SED), the appearance of disks, temperature profile, and chemistry. 
Dust coagulation is also the first step towards planet formation. 

Aims. We investigate dust growth due to settling in a ID vertical column of a disk. We gradually build up the complexity 
of the models by considering the effects of porosity, different collision models, turbulence, and different gas models 
respectively. This way we can distinguish the effects of these physical processes on particle growth and motion. It is 
known from the 10 micron feature in disk SEDs, that small micron-sized grains are present at the disk atmosphere 
throughout the lifetime of the disk. We hope to explain such questions as what process can keep the disk atmospheres 
dusty for the lifetime of the disk and how does the particle properties change as a function of height above the 
midplane. 

Methods. We use a Monte Carlo code to follow the mass and porosity evolution of the particles in time. The used 
collision model is based on laboratory experiments performed on dust aggregates. As the experiments cannot cover all 
possible collision scenarios, the largest uncertainty of our model is the necessary extrapolations we had to perform. We 
simultaneously solve for the particle growth and motion. Particles can move vertically due to settling and turbulent 
mixing. We assume that the vertical profile of the gas density is fixed in time and only the solid component evolves. 
Results. We find that the used collision model strongly influences the masses and sizes of the particles. The laboratory 
experiment based collision model greatly reduces the particle sizes compared to models that assume sticking at all 
collision velocities. We find that a turbulence parameter of a = 10~^ is needed to keep the dust atmospheres dusty, 
but such strong turbulence can produce only small particles at the midplane which is not favorable for planetesimal 
formation models. We also see that the particles are larger at the midplane and smaller at the upper layers of the disk. 
At 3-4 pressure scale heights micron-sized particles are produced. These particle sizes are needed to explain the 10 
micron feature of disk SEDs. Turbulence may therefore help to keep small dust particles in the disk atmosphere. 

Key words. Methods: Numerical, Planets and satellites: formation, Protoplanetary disks 



1. Introduction 

The coagulation of dust particles is the first step towards 
planet f ormation. As dust par ticles are the main source of 
opacity (jSemenov et al.l . I2003D , they also determine many 
observable quantities of disks, such as the spectral energy 
distribution and scattered light images (see below). Due to 
the vertical component of the stellar gravity, dust particles 
sediment towards the midplane of the disk. The relative set- 
tling velocity between particles of different aerodynamical 
properties drives coagula tion and this process was investi- 
gated by several authors (iSafronovl 1969t Nakagawa et al 



1986 : ISchrapler fc Hennina. 



2004 iDullemond k Dominikl . l2005[ ). 



2004 Dullemond fc Dominik 



Observational evidence of the vertical sedimentation of 
grains exists for a large nu mber of disks, although su ch ev- 
idence is usually indirect (jHenning fc Meeua 1201 it) . Sub- 
micron grains are present at the disk surfaces as shown 
by scattered light images in the optical and near infrared 
(NIR) wavelengths. Such multi-wavelength scattered light 



image s provide evidence for g rain growth (IWatson et al.l 
(|2007D and references therein). iPinte et aD (|2007[ 1 showed 
by reproducing multi-band images of the binary system of 
GGTau that the dust scale height for 10 micron-sized par- 
ticles is roughly half of that for micron-sized particles. The 
spectral energy distribution (SED) is also affected by set- 
tling. D'Alessio et al.l ()2006t ) showed that in order to ex- 
plain the median SEDs of classical TTauri stars, the dust 
to gas ratio has to be reduced by a factor of 10 at the 
disk atmosphere compared to the standard value. There are 
also indications that the settling of grains i s corr elated with 
the age of the disk (jSicilia-Aguilar et al.l . l2007l ). However, 
the connection between the exact shape of the 10 mi- 
cron feature of SEDs and sedimentation is not well under - 
stood dP uUcmond & Dominik, 2008; Juhas z et all I2010D . 
lApai et al.l (|2005t ) showed evidence for settling in disks 
around brown dwarf stars. They concluded that growth, 
crystallization and settling of dust happens around low 
mass stars in a similar manner as around intermediate and 
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solar mass stars. This suggests that the first stage of planet 
formation is a robust process occurring in all kinds of disks. 

Sedimentation also affects the ve rtical temperature 
struct ure of the disk. The simulations of lAikawa fc Nomural 
(|2006(l showed that as the dust particles sediment towards 
the m.idplane, the opacity is reduced, and the temperature 
of the gas decreases. As the stellar radiation can now pen- 
etrate deeper in the disk, the temperature at intermediate 
heights increases. The change in the density and tempera- 
ture structure as well as the radiation field n aturally influ- 
ences the chemistry of the disk atmospheres (jBergin et al.l . 
[2OOI . 

Recently IVasvunin et al.l ()2011h investigated the effects 
of dust growth on disk chemistry and found that the main 
effect of coagulation on the disk chemical composition 
comes mostly from sedimentation in their models. 

Dust sedimentation not only affects the upper layers 
of the disk, but also the midplane regions. As long as the 
dust to gas ratio is much smaller than unity, the dust can 
be treated as a passive tracer in the gas flow. However, if 
the dust to gas ratio becomes comparable to unity, the dust 
will influence the gas. The dust can be accumulated by sed- 
imentation around the midplane of the disk. In such a sce- 
nario, a shear is present between the dusty midplane layer 
and the less dusty upper layers. This shear triggers Kelvin- 
Helmholtz instability which de velops into turbulence. This 
process was first recognized bv I Weidenschilline (Il980l) and 
is stil l under a c tive i n vestigation (e.g. by iJohansen et al.l 
(1200611 ■ IChianel (l2008l l. ITee et al.l (120100 ). 

A collaboration started between the lab-community and 
the modelers to better constrain the dust evolution in pro- 
toplanetary disks, using a realistic collisio n model that 
is bas ed on the laboratory experiments. In iGiittler et al.l 
(|2010D (h enceforth Paper I ), we introduced this collision 
model. In lZsom et al.l ()2010[ ) (henceforth Paper II), we used 
this collision mod el for the first time i n the Monte Carlo 
(MC) method of IZsom fc Dullemondl (|2008D (henceforth 
ZsD08). The models of Paper II were local box models, 
meaning that the dust evolution was only followed at one 
location of the disk. These models showed that bouncing 
plays an important role in dust evolution. 

We further develop these models to simulate a ID ver- 
tical column in the disk, thus investigating sedimentation- 
driven coagulation. We want to better understand the pro- 
cess of sedimentation and the role of particular physical 
phenomena like porosity of the aggregates, collision models 
and turbulence. 

Previous work bv lDullemond fc Dom inik (2005) (hence- 
forth DD05) showed that without a mechanism that re- 
duces the sticking probability of particles in the upper lay- 
ers of the disk, or without a continuous source of small 
particles, the observed spectral energy distributions (SED) 
of TTauri stars should exhibit a very weak IR excess. In 
contrast, the ob served SEDs of TTauri stars have strong 
IR ex ce sses (e.g. Furlan et al.l ( 20051 ) ; iKessler-Silacci et all 



(|2006() : iBouwman etal 



()2008 )). Therefore some gram- 
retention mechanism is needed to explain the SEDs. 

Previous models of grain evolution assumed a contin- 
uous cycle of growth and fragmentation, w hich provides 
the ne c essary amount o f sma ll particles (e.g. iBrauer et al.l 
(|2008[) . iBirnstieletaP (|2009D V However, Paper I and II 
showed that particle evolution is halted by bouncing and 
no cycle of growth and fragmentation is present. In this 
paper we simulate dust evolution driven by Brownian mo- 



tion, turbulence, and sedimentation in a ID vertical col- 
umn of the inner disk and the additional physics of Paper 
I and II are included. We investigate the time evolution 
of sedimentation-driven coagulation, and search for ways 
that can keep a sufficient amount of the small dust par- 
ticles levitated at several pressure scale heights to explain 
the observed SEDs of young stars. 

The paper is organized as follows. In Sect. [5] we de- 
scribe the numerical method used to follow the particle 
motion and coagulation. We validate the code and increase 
the complexity of the model step-by-step in Sect. [31 We 
show the results in Sect. JH finally we discuss those results 
in Sect. [5] and provide a summary in Sect. [HI 



2. Numerical method 

2.1. Basic considerations 

The local box approach (or "particle-in-a-box" approach) 
in Paper II is based on two assumptions. 1.) The particles 
are homogeneously mixed inside the box. 2.) Particles do 
not enter or leave the box, i.e. it is closed. Due to these 
two assumptions, it was not necessary to follow the exact 
location of the particles. 

In the models considered here, however, we place such 
boxes (or grid cells) on top of each other to simulate a ID 
column in the disk and follow how particles settle towards 
the midplane. Inevitably, particles move from box to box 
during this process. Therefore the assumption that particles 
cannot enter or leave the boxes has to be relaxed. The first 
assumption of the method in Paper II is kept, we still as- 
sume that during the coagulation calculation, the particles 
inside a given box are homogeneously distributed (however, 
for particle motions, the individual positions of the parti- 
cles are used). The second assumption is modified in the 
following way. 2.) The simulated column is closed, e.g. par- 
ticles inside the column can move freely vertically. However 
neither do new particles enter from the "outside" , nor do 
particles from inside the column leave. As particles move 
through the boxes, it is necessary to follow the position of 
the particles (see Sect. 12. 3|) as we must find out in which 
box a particle is located. 

The motion of particles imposes a limit on the time step 
of the simulation. We do not want the particles to move 
more than one box in a time step. A sedimenting particle 
should have the possibility to interact with all other parti- 
cles along its way, it should not skip over boxes thus avoid- 
ing the particles in it. Therefore, we use an adaptive time 
stepping method. The maximum of all particle velocities is 
obtained (wmax), and since we know the height of the boxes 
(^box), the maximal (safe) time step can be determined as 



At^C- 



(1) 



where C is the Courant number which we typically set to 
be 0.1. 

The code schematically performs the following steps: 

1. The velocities of the particles are calculated. 

2. A safe time step is determined to avoid particle 'jumps'. 

3. The position of the particles is updated using their ve- 
locities, their previous positions, and the time step. 

4. We determine the box in which each particle resides. 
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5. We call the coagulation subroutine described in Paper II 
to calculate the evolution of the particles separately in 
each box for the given time step. There can be multiple 
collisions per time step. 

2.2. Initial conditions 

We assume that the gas density profile is constant in time 
during the simulation. This assumption is valid if the sim- 
ulated time is less or comparable to the viscous timescale 
of the gas. The viscous timescale can be calculated as 



iv 



r^/vT, 



(2) 



where r is the distance from the central star, vt is the tur- 
bulent viscosity. The value for tvis at 1 AU varies between 
10'^ and 10^ yr depending o n vt- We assume that turb u- 
lence is parameterized by the IShakura fc SunvaevI (|1973[ ) a 
parameter 



i^T = aCsH„ 



where Cg is the isothermal sound speed, and Hg is the pres- 
sure scale height of the gas disk. The turbulence parameter 
a reflects the strength of the turbulence in the disk. Typical 
values range between a — 10^® and 10^^, where the for- 
mer corresponds to the turbulent strength in dead zones, 
the latter describes turbulence in disk atmospheres. In this 
paper, we assume that a is constant as a function of height, 
which may be changed in future work (see Sect. 15. 3p . 

The vertical structure of the disk is determined by the 
equilibrium between the vertical component of the gravita- 
tional force and the vertical pressure gradient in the gas. 
If the disk mass (Mdisk) is much smaller than the mass of 
the star (M*), and the vertical thickness of the disk (Hg) is 
a small fraction of the radial distance (both conditions are 
safely met for the disk parameters described below), then 
the vertical density can be approximated as 



Pgir,z) 



S(r) 



2TrH, 



exp(-zV2ij2 



(4) 



where S(r) is the gas surface density at distance r, and z 
is the height above the midplane. In this paper we choose 
M* = O.SMq, r = 1 AU, i;(l AU) = 100 g/cm^ similarly 
to DD05. The pressure scale height can be calculated as 



Hg = Cs/Q, 



(5) 



where il is the orbital frequency at 1 AU. The isothermal 
sound speed is 






(6) 



where fc^ is the Boltzmann constant, fi is the molecular 
weight, which is 2.3 for molecular gas, mp is the mass of 
the proton, and T is the temperature of the gas, which is 
200 K for the stellar and disk parameters considered above 
(see DD05). We assume the temperature to be constant 
as a function of height. This is a reasonable assumption 
well below the photosphere if the temperature of the gas is 
solely determined by the stellar irradiation and the thermal 
coupling to the grains. The height of the photosphere may 
(and presumably will) change as the particles rain out, but 
we do not include this effect in this paper. 

We need to quantify the proper number of particles and 
boxes to be used. As we know from ZsD08, a low number of 



particles is not desirable because the physics of the collision 
kernel will not be reproduced properly. However, many par- 
ticles make the code computationally expensive and unfea- 
sible to run. If too few boxes are used, the Gaussian profile 
of the gas disk will not be resolved properly. If too many 
boxes are used, sufficient number of particles per box is 
needed, which again renders the simulation computation- 
ally expensive. We found by numerical experiments that 
using 10^ particles and 40 evenly spaced boxes is a good 
compromise. We simulate 4 pressure scale heights (0.16 AU 
above the midplane), therefore there are 10 boxes per pres- 
sure scale height to resolve the Gaussain profile of the gas. 
10^ particles are also enough to remove numerical artifacts 
from the simulation and reproduce the physics of the colli- 
sion kernel properly. Using 10^ particles and 40 boxes also 
mean that there is at least 1 particle in the uppermost box 
in the beginning of the simulation. 



(3) 2.3. Position update 



The vertical position of the particles are determined by ver- 
tical settling and turbulent diffusion. In principle, Brownian 
motion also contributes to the change of particle positions, 
but its effect is negligible compared to the other two effects. 
The equation governing the diffusion and settling 
of the dus t in a non- h omogenous gas densit y field 
is (DubruUelian, 119951: iFromang fc Papaloizoul 120061: 
iCieslai . i201Ql 



dpd 
dt 



d_ 

dz 



DdPg 



d 



+ -^{nhspdz), (7) 

oz 



or equivalently: 



dpd d f dpd\ d f 1 dpg\ d 2, n 

-m=d-z r^^j-a; r " '"'y.^rd-z^''''''' ''^ 

(8) 
where pd is the dust density, Dd is the diffusion coefficient 
of the dust (see the next paragraph) and ts is the stopping 
time of the particle. The stopping time is the timescale a 
particle needs to react to the changes of the surrounding 
gas. We define the dimensionless Knudsen number as 



Kn = 



X 



nifp 



(9) 



where a is the size of the aggregate, and A,nfp is the mean 
free path of the gas. A particle is in the Epstein regime 
if Kn < 1 (to b e more p recise, if a < |Amfp), where the 
stopping time is (JEpsteinl (^1924.)): 



ts — ^Ep — 



3m 



iVthPgA' 



(10) 



where m and A are the mass and the aerodynamical cross- 
section of the particle, and uth is the thermal velocity. If the 
Knudsen number is greater than 1 (at high gas densities 
where the mean free path is low or in the case of large 
particles), the Stokes regime applies and the stopping time 

becomes 

3m A a ,, , 

t, = tst^ tX . (11) 

AvthPgA 9 Amfp 

The first term on the right hand side of Eq. [8] is the 
diffusion term. Using only this term, particles with tg = 
(tracers) would be homogeneously distributed as a function 
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of height over several diffusion timescales. The first and 
the second term together on the right hand side ensures 
that the tracer particles will be distributed according to 
the background gas density field. The third term describes 
the settling of the particles. Equation[5]is valid if the motion 
of the dust does not influence the motion of the gas (the 
back-reaction from the dust to the gas is negligible). This 
condition is met if the dust to gas ratio is ^ 1. 

We note that we do not solve for the dust density di- 
rectly. We follow the motion of dust particles, each of which 
represents a portion of the total dust mass inside the col- 
umn. Thus we derive the corresponding velocities (or fluxes) 
for the first, second, and third terms of Eq. [5]to calculate 
the position update of the particles. 

We calculate Dd, the diffusion coefficient of the dust, 
and define the diffusion velocity, f pi . The diffusion coef- 
ficien t of the gas can be defined as (jShakura &: Sunvaevl . 
[19731) 



D„ 



Ut 



aCsHg. 



(12) 

Based on lYoudin fc LithwicM ()2007D . the diffusion coeffi- 
cient of the dust can be calculated as 



where St is the Stokes number 

St = tM. 



(13) 



(14) 



The average displacement of a particle in ID during the 
time step of At then is 



L = y/2DdAt. 



(15) 



The real displacement of the particle (Az) is drawn from 
a Gaussian distribution which has zero mean and a half 
width of L. The "diffusion velocity" can then be calculated 
as 

VDi = Az/At. (16) 

This velocity component tries to smear out dust concentra- 
tions. It is important to note that the real, physical velocity 
of the particle during turbulent diffusion changes randomly 
every time the aggregate interacts with a turbulent eddy. 
The diffusion velocity defined above is a numerical con- 
struct to calculate the time-averaged velocity of the particle 
during a time-interval At. 

The second term on the right side of Eq. |S] results in a 
systematic velocity term which pushes particles towards the 
density maxima of the gas. This velocity can be determined 

by 

VD2=Dd-^. (17) 

Pg dZ 

Using these two velocity components (wdi and vd2), the 
particles will be distributed according to the gas density 
profile in a diffusion timescale. 

The fact that particles with tg > settle towards the 
midplane and have a scale height less than Hg is the result 
of the third term of Eq. [8l the settling velocity. The settling 
velocity of a particle is given by 

t>sct = -zQhs. (18) 

The new position of the particles can then be deter- 
mined by using these three velocity terms: 

Z = 2old + {VDI + VD2 + Vsct)At. (19) 

This description os identical to the one used in ICieslal 

(|2nifiD . 
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Fig. 1. The collision rate between dust aggregates as a func- 
tion of particle mass to particle mass measured one simula- 
tion of Paper II (the MMSN model with ID Mtld-4ml00). 
Initially particles grow due to Brownian motion and col- 
lisions between equal sized aggregates are frequent. When 
turbulent relative velocity dominates over Brownian motion 
(masses above 10* g), collisions between equal sized aggre- 
gates are less frequent and collisions between different sized 
aggregates are dominant. 



2.4. Coagulation 

The collision model used in this work is similar to the one 
used in Paper II. There are however two differences. The 
first difference is the additional source of relative velocity 
due to differential vertical settling (see Eq. [T5| : 



Avs = ll^sotl - Wsct2|- 



(20) 



The second difference concerns the calculation of the 
aerodynamical cross section which is used to calculate the 
stopping time. In Pape r II we used the ge ometrical cross 
section of the particles (jOrmel et al.L l2007f) 



A 



2^VI/2/3, 



(21) 



where r^ is the compact radius of the aggregate (assuming 
that the mass of the particle is contained in a compact 
sphere of radius Tc), and ^ is the enlargement parameter. 
\1/ = 1 for co mpact particles and the value is higher for 
fluffy particles (jOrmel et al.L 120071 ). This formula works well 
for particles with fractal dimensio n above 2. H owever, we 
use the porosity model of Okuzum i et al.l (|2009" ). and their 
model produces aggregates with fractal dimension below 2. 
If one calculates the stopping time of such an aggregate in 
the Epstein regime (Eq. [T0| using the formula in Eq. [211 
one gets that the stopping time is less than the stopping 
time of a monomer. This is clearly unphysical. The reason 
for this low stopping time (large area) is that Eq. [21] does 
not take into account the empty space between the 'fractal 
branches'. To avoid such unphysical results for aggregates 
with low fractal dimensions, we use the aerodyn amical cross 
section as defined in Eq. 47 in lOkuzumi et al.l (|2009iV 

2.5. Porosity and collision models 

It is generally assumed that dust particles in the proto- 
planetary disk are aggregates, i.e. they consist of micron- 
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Table 1. Overview and results of all the sedimentation simulations. 



ID 


[g/cm^] 


Coll. model 


Por. model 


a 


*max 


< 


*fln > 


[g] 


< St> 




(1) 


(2) 


(3) 


(4) 


(5) 


(6) 




(7) 


(8) 


(9) 


(10) 


DD 


100 


hit&stick only 


compact 





1 




1 


2.5 X 10-" 


3.5 X lO-' 





DDa 


100 


hit&stick only 


Okuzumi 





10^ 


1.1 


_ X 10^ 


1.8 X 10** 


1.6 





SBl 


100 


sinipl. Br. 


Okuzumi 





10* 


2 


X 10^ 


2.2 X 10"=^ 


8.0 X 10-^ 





SB2 


100 


simpl. Br. 


Okuzumi 


io-« 


10* 


i.e 


) X 10^ 


3.2 X 10-^ 


7.0 X 10-^ 


0.2 


SB3 


100 


simpl. Br. 


Okuzumi 


10-* 


5x 10^ 




21 


2.2 X 10-* 


5.0 X 10-* 


0.6 


CBl 


100 


compl. Br. 


Okuzumi 


10-* 


10^ 




3 


1.2 X 10-" 


4.5 X 10-^ 


0.2 


CB2 


100 


compl. Br. 


Ormel 


10-* 


10 




2 


1.2 X 10-" 


4.7 X 10-^ 


0.2 


CB3 


100 


compl. Br. 


Okuzumi 


10-" 


3x 10" 




3 


3.6 X 10-^ 


1.6 X 10-* 


0.95 


CB4 


100 


compl. Br. 


Ormel 


10-" 


10 




2 


8.9 X 10-* 


1.2 X 10-* 


0.95 


LDl 


10 


compl. Br. 


Okuzumi 


10-" 


5x 10" 




2 


8.9 X 10-" 


9.7 X 10-^ 


0.2 


LD2 


10 


compl. Br. 


Okuzumi 


10-" 


5x 10' 




23 


2.3 X 10"^ 


1.3 X 10-* 


0.95 


HDl 


1000 


compl. Br. 


Okuzumi 


10-* 


5x 10* 




6 


8.3 X 10-" 


6.7 X 10-* 


0.5 


HD2 


1000 


compl. Br. 


Okuzumi 


10-" 


10=^ 




3 


1.2 X 10-^ 


4.6 X 10-^ 


0.95 



Col. 1 is the ID of the simulations (DD - Dullemond& Dominik models, SB - simplified Braunschweig model, CB - complete 
Braunschweig model, LD - low density model, HD - high density model) , col. 2 is the gas surface density, col. 3 describes the used 
collision model ( hitfcstick, simplifie d B raunschweig model, or complete Braunschweig model), col. 4 indicates the used porosity 
model (based on lOrmel et al.l (|2007D or lOkuzumi et aLl (|2009D ). col. 5 describes the value of the turbulence parameter a, col. 6 is 
the maximum enlargement parameter of the aggregates during the simulation, col. 7, 8, and 8 are the final average enlargement 
parameter, mass and Stokes number of the particles respectively, and col. 10 is the scale height of the dust expressed in the scale 
height of the gas. 



sized sphe res (monomers) which are co nnected via sur- 
face forces iHenning fc Stoenienkol (|1996[ ). The material of 
these monomers can be iron, silicate, organics (tar), differ- 
ent ices and a mixture of these (e.g., silicates with organic 
or icy mantels). The collision model described in Paper I 
is based on laboratory experiments performed on room- 
temperature dust aggregates made out of micron-sized sili- 
cate monomers. This naturally limits the region in the disk 
where our collision model is applicable. Therefore, our re- 
sults are not applicable in regions where ices condensate 
(beyond the snow-line) or close to the central star where 
the effects of sinteri ng (partial me lting and crystallization) 
becomes important ()Poppel . l2Q03D . Another limitation that 
originates from the adopted collision model is that the labo- 
ratory experiments were performed using fractal dimension 
3 aggregates. Therefore, we cannot reliably consider the ini- 
tial fractal growth of the aggregates. This topic is further 
discussed in Sect. 15.11 

The monomers a re embedded in t he ga s and move rel- 
ative to each other iBeckwith et al.l (|200Cll ). The particles 
collide, stick, and grow due to their relative motions. There 
are two limiting cases for growth: cluster-cluster aggrega- 
tion (CCA) and particle-cluster aggregation (PC A). CCA 
growth means that the aggregates tend to collide with other 
aggregates that consist of the same number of monomers. 
PCA growth occurs when an aggregate collides successively 
with much smaller aggregates (monomers) only. The re- 
sulting structure of the aggregates are different for the two 
growth types: CCA growth produces fluffy fractal aggre- 
gates with a fractal dimension of 2, while PCA growth 
pr oduces compact structures with a fractal dimension of 
3 (lOssenkonJ . [19931) . 

Two monomers can perform different types of motions 
while they are connected if the available kinetic energy is 
high enough. These are rolling, slidin g, twisting, and the 
connection can break due to pull-off (iDominik fc Tielensl . 



'1997*). Most important of these is rolling as the aggregate 
can dissipate the initial collisional energy efficiently through 
rolling motions. The rolling energy (i?rou) is defined as the 
energy that is needed to roll two monomers by 90 degrees. 
During the initial stages of growth, the collisional energy 
is lower than the rolling energy, therefore no significant re- 
structuring occurs. The aggregates stick at the first contact 
and we refer to this as the hit&stick (SI) collision type. 

We adopt the hit&stick porosity model of 
lOkuzumi et al.l (|2009t) in thi s paper. In Paper II we 
used the porosity model of lOrmel et al.l (|2007D which 
only treats PCA and CCA collisions and constructs 
semi- analytical recipes if the size (or mass) ratio of 
the two collidi n g agg regates are intermediate. However, 
lOkuzumi et al.l ()2009() improved on this by modeling 
these intermediate regions, so called quasi-CCA collisions 
(QCCA), where two clusters with a given mass ratio 
can collide. The porosity model of Ormel et al. produces 
aggregates with a fractal dimension of 2.5. The Okuzumi et 
al. model results in fluffy, fractal dimension 2 aggregates. 
As the Ormel model analytically interpolates between 
PCA and CCA collisions, whereas the Okuzumi et al. 
model directly simulates these collisions, we prefer the 
Okuzumi model as the fiducial hit&stick porosity model. 

As the particles grow, their collisional energy is in- 
creasing and we cannot neglect restructuring. This phase 
starts when the collisional energy is a few times larger than 
the rolling energy between the monomers. As compaction 
occurs, the fractal dimension of the aggregates increases. 
At the same time, the average number of connections a 
monomer has (coordination number) increases as well. 

The further evolution of the dust aggregates is as of yet 
unclear and a matter of ongoing debate. Currently there 
are two competing collision models. The Braunschweig col- 
lision model (Paper I) is mostly based on laboratory exper- 
iments performed by using fractal dimension 3 aggregates 
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(so called dust cakes and their pieces) . As such, it postulates 
that the aggregates at one point during their evolution in 
the protoplanetary disk reach a rather compact state with a 
fractal dimension of 3. The validity of this assumption was 
not thoroughly checked yet (but see future work in Sect. 
15. ip . This collision model predicts that nine different col- 
lision types can occur between dust aggregates of different 
sizes and porosities (e.g., bouncing, mass transfer, crater- 
ing, etc.). 

The other alternative c omes from molecular dynamics 
models. 'Wada et al. (2008) performed a series of collisions 
using equal sized aggregates of fractal dimension 2 with an 
ever increasing relative velocity. They found that the aggre- 
gates end up with fractal dimension 2.5 due to restructuring 
before fragmentation sets in. This would imply that the col- 
lision model can be constructed from hit&stick collisions at 
low velocities, sticking with restructuring at higher veloc- 
ities, and finally fragmentation. They did not observe any 
other collision type. 

High relative velocities are required to initiate restruc- 
turing of equal sized collisions. Relative velocity due to 
Brownian motion is negligible for large aggregates, as Avb 



decreases with mass. Relative velocity due to differential 
radial drift and settling is zero, as equal sized aggregates 
drift and settle with the same speed. Relative velocity due 
to turbulence is also zero or negli gibly small, if both parti- 
cles are well coupled t o the gas (jWeidenschilling fc Cuzzi 
Il993l : lOrmel fc Cuzzil . l2007f) . If one particle is decoupled 
from the gas, the relative velocity between different and 
equal sized aggregates is roughly constant. However, one 
needs sufficiently large aggregates to enter this regime. 

The relative velocity between different sized aggregates 
is larger than between equal sized aggregates during the 
initial stages of growth. As a result, collisions between dif- 
ferent sized collisions are more frequent. To support this 
statement, we use one simulation result obtained in Paper 
II: one of the Minimum Mass Solar Nebula (MMSN) mod- 
els with ID Mtld-4ml00. We produced a logarithmically 
spaced mass grid array and we calculated how many col- 
lisions happened between the aggregates within each grid 
cell during the simulation. The number of collisions per 
grid cell is indicated by the contour levels in Fig. [TJ We 
only used the first 300 yr of the simulation because par- 
ticles at i > 300 yr experience other collision types than 
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Fig. 3. Mass distribution for the DDa model using the Okuzumi hit&stick porosity modeh The axes and contours are 
the same as in Fig. [51 



sticking. Since the relative velocity of small particles is 
dominated by Brownian motion, similar sized aggregates 
collide (CCA-like growth) as long as the particle mass is 
less than 10~^ g. However, when the turbulent relative ve- 
locity dominates over Brownian motion (at masses above 
10~^ g), the collision rate between equal sized aggregates 
is reduced. The large aggregates (mass of 10""' g) preferen- 
tially collide with smaller ones (mass of 10~* g). Therefore 
the results of Wada ct al. ( 2008) (which is based on equal 
sized collisions) might not capture the full complexity of 
the aggregate restructuring phase. 

Ideally, a restructuring collision model considering dif- 
ferent sized aggregates is needed. As such a model - and 
its experimental verification - is not available yet, we for 
simplicity decide to stick with the Braunschweig lab-based 
model, as we did in previous works. We acknowledge the 
potential caveats of our adopted collision model, which we 
further discuss and quantify in Sect. 15.11 

3. Prelude: the effects of porosity and collision 
models 

We perform altogether 21 simulations to investigate the ef- 
fects of different porosity models and turbulence, in which 



we gradually use more realistic collision models. The IDs 
and parameters of these simulations are shown in Tab. [TJ 
First we compare our model against the results of DD05 for 
compact particles (model DP in Tab. [H. Th e n we use the 
hit&stick porosity model of lOkuzumi et all (|2009[ ) to in- 
vestigate the effects of porosity (model DDa). So far we as- 
sume that the aggregates stick at all relative velocities and 
that the turbulence parameter a is zero. In the next step 
we construct a more realistic collision model with sticking, 
bouncing and fragmentation (model SBl). We call this col- 
lision model the "simplified Braunschweig model" because 
it uses only three collision types out of 9, which is described 
in Paper I (the complete Braunschweig model). In the next 
step we turn on turbulence (models SB2-3) to examine the 
effects of turbulent stirring. Finally, we use the complete 
Braunschweig model with turbulence (models CBl-4) and 
use different disk models (models LDl-2 and HDl-2). 

3.1. Test: comparison with the DD05 model 

DD05 performed a simulation (S2 in their paper, DD in 
this paper), where the disk model is the same as the one 
described in Sect. 12.21 the particles are compact, upon col- 
lision particles stick together at all collision energies, and 
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the only source of relative velocity is Brownian motion and 
differential settling. They found that the 'rain-out' parti- 
cles reach the midplane in 500 yr having attained sizes of 
a millimeter (10~^ g in mass). 

We performed the exact same simulation to validate our 
code. We find that the 'rain-out' particles in our simulation 
reach the midplane also in 500 yr and have masses of 10^^ 
g. Therefore we can conclude that our code works properly. 
We illustrate the mass evolution of particles as a function 
of their height at six different snapshots {t — 1 yr, 10 yr, 
100 yr, 316 yr, 10^ yr, 10^ yr) in Fig.[ll 

Brownian motion is essential in our simulations because 
growth by Brownian motion initializes the sedimentation- 
driven coagulation. The reason is that we have initially a 
mono-disperse particle size-distribution (meaning that all 
particles have the same size and mass). The aerodynami- 
cal properties of these particles are identical, thus there is 
no relative velocity due to settling between the monomers 
at a given height. If growth due to Brownian motion was 
not initiated (e.g. growth by Brownian motion did not in- 
troduce aggregates with different aerodynamical properties 
than that of the monomers), the monomers would simply 
sediment to the midplane without any growth. Although 



DD05 included Brownian motion, this effect would not have 
been present as that simulation started with a (narrow, but 
not infinitely narrow) size distribution. 



As shown in Fig. [21 growth by Brownian motion is faster 
at the midplane due to the higher gas and dust densities (at 
t = 1 and 10 yr). Once particles in the upper layer also start 
to grow by Brownian motion, sedimentation-driven coagu- 
lation starts and particles at the upper layers grow much 
faster than the aggregates at the midplane (at t = 100 
yr) because the absolute value of the settling velocity in- 
creases with height (see Eq. [TH]). The heaviest particles 
sweep up the smaller particles while they sediment and fur- 
ther increase their settling velocity, resulting in a rain-out 
at i = 500 yr. Once the first rain-out particles reach the 
midplane, they could only grow by Brownian motion, be- 
cause at the midplane, the settling velocity of any particle is 
zero (see Eq. [T51) . But the relative velocity due to Brownian 
motion for such heavy aggregates is negligible. Therefore, 
the particles that have reached the midplane, do not in- 
crease in mass any longer. 
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3.2. The effects of the hit&stick porosity model 

In the previous section we used compact particles. In this 
section, we consider dust aggregates which are built from 
(sub)micron-sized solid spheres, called monomers. Such 
monomers collide with each other and form fluffy struc- 
tures due to hit-and-stick collisions (collisions that result in 
sticking upon the first contact). In this section we assume 
that all collisions result in sticking. We include growth by 
Brownian motion and settling only. 

We use here th e hitfc stick porosity model constructed 
bv lOkuzumi et al.l (|2009[ ). As discussed in Sect. [131 this 
collision model defines a third type of aggregation next to 
the PCA and CCA collisions, that is QCCA (quasi cluster- 
cluster aggregation - collisions between clusters with a pre- 
defined mass ratio). This model is based on simulations of 
hit&stick collisions without restructuring. 

The mass and the porosity evolution are shown in Figs. 
[3] and m The most striking property of this simulation is 
the maximum mass and porosity of the particles. We end up 
with particles of 10^" g in mass having an enlargement pa- 
rameter of almost 10* (while the compact radius of such an 
aggregate is some meters, the enlarged radius is several kilo- 



meters!). The stopping time in the Epstein regime fEa.[TU| 
is proportional to m/A. As the Okuzumi-model produces 
aggregates with a fractal dimension of ~ 2 (the mass scales 
with a^, where a is the particle radius), the stopping time 
only slightly increases with mass. Therefore, particles set- 
tle slowly and produce extremely fluffy structures. At one 
point, however, the aerodynamical cross section radius be- 
comes larger than the mean free path of the gas, and the ag- 
gregate enters the Stokes regime (Eq. [TT|) . As the stopping 
time is now proportional to ma/ A, the stopping time more 
strongly increases with mass. This transition from Epstein 
to Stokes regimes happens at t=900 yr for the particles lo- 
cated at 1.7 Hg above the midplane. Once the transition 
happens for a given particle, it settles to the midplane in 
a matter of years due to the heavy mass of the aggregate. 
Therefore, when the size of the aggregate becomes compa- 
rable to the mean free path of the gas (a = Amfp), we reach 
a natural upper size limit where rain-out in any model is 
expected. 

We emphasize that any collision model containing ex- 
clusively sticking is only valid, if no significant restructur- 
ing happens during the collisions (e.g. the collision energy 
is less than 5 £"1011, the rolling energy, which is the energy 
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needed to roll two monomers by 90 deg). This condition is 
clearly not met at all times in our simulations, e.g. the rain- 
out particles can have collision velocities with the swept up 
particles as high as several 10 m/s in these simulations. 
Such collisions would result in catastrophic fragmentation. 
Therefore, the results presented in this section should be 
considered as toy models. 

We conclude here that porosity alone is not sufficient to 
explain the observations. Therefore, the hit-and-stick as- 
sumption is insufficient to describe all settling stages, since 
we expect that the condition E < 5i?roii will be broken at 
some point. Clearly, a more realistic collision model that in- 
clude compaction and fragmentation of dust aggregates, is 
then needed. Next, we will investigate the consequences of 
including these physical regimes, first by using a simplified 
prescription. 

3.3. A simplified Braunschweig model 

In this section we use the simplified version of the collision 
model described in Paper I. We assume sticking and the 
increase of the porosity, if the collision energy is smaller 
than 5i?ron- Bouncing with compaction is used if the colli- 
sion energy is greater than SE'roii, but the relative velocity 



of the two aggregates is less than 1 m/s. Fragmentation 
occurs if the relative velocity of two aggregates is greater 
than 1 m/s. The recipe for mass and porosity evolution 
for bouncing and fragmentation is taken from Paper I (our 
hit & stick (SI), bouncing with compaction (Bl), and frag- 
mentation (Fl) collision types). We still assume that the 
particles grow by Brownian motion and settling only (the 
effects of turbulence is discussed in the next section), and 
we use the porosity model of Okuzu mi et al. (2009) for the 
hit&stick phase. 

The evolution of the mass is shown in Fig. [51 The mass 
distributions a.t t = 1, 10, 100 yr are identical to Fig.|31 At 
t = 316 and 1000 yr we see the effects of bouncing at the in- 
termediate energies. The rain-out particles cannot increase 
their mass further, when they suffer bouncing collisions. 
Therefore, their masses are 10~^ g when they arrive to the 
midplane and no larger aggregates are formed as opposed 
to the DDa model. As the rain-out particles are smaller, 
they settle slower, and reach the midplane only at i = 6000 
yr. The enlargement parameter at the end of the hit&stick 
phase is ^ 10^ and it decreases significantly only in the 
midplane to values between lO-lO** as the rain-out particles 
collide and bounce with the aggregates in the midplane. 
The particles in this simulation are always in the Epstein 
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regime, as bouncing collisions prevent growth to sizes above 
the mean free path of the gas. 

We conclude that realistic collision models reduce the 
particle sizes in sedimentation-driven coagulation and thus 
reduce the settling time. However these collision models 
with the combined effects of porosity are still insufficient to 
explain the long-term presence of the 10 micron feature in 
disk SEDs. We now include turbulence into the simulation. 



4. Results 

4.1. The efFects of turbulence 

So far all particles sooner or later ended up at the midplane 
because there was no effect that could counteract settling. 
In this section we examine the effects of a non-zero turbu- 
lence parameter, which can stir particles back up. 

A small turbulence parameter (a = 10~^) does not sig- 
nificantly affect the masses of the rain-out particles com- 
pared to the SBl model of the previous section (see model 
SB2 in Tab. [T|). As particles do not only settle but also 
diffuse downward (and upward) due to turbulence, the 
time the rain-out particles reach the midplane is some- 
what shorter than for model SBl. In all previous simula- 



tions a dense dust layer formed at the midplane of the disk. 
However, even this low level of turbulence can prevent the 
formation of this layer and introduce a non-zero (although 
small) dust scale-height. The porosity of the aggregates are 
also similar to the results described for the SBl simulation. 

The infiuence of turbulence is more pronounced if a = 
10^^. The mass evolution of the aggregates is shown in 
Fig. [S] The first rain-out particles reach the midplane al- 
ready at i = 500 yr due to downward diffusion, although 
these particles have lower masses than in model SBl (10~^ 
g - therefore, in the absence of turbulence, these particles 
would reach the midplane later than the particles in SBl). 
The porosity of the aggregates is strongly influenced by the 
higher turbulence value and bouncing. Due to the increased 
turbulent relative velocities, particles start bouncing before 
they reach the midplane. The average enlargement param- 
eter is 2 X 10'^ at the end of the hit&stick phase. As the 
particles settle to the midplane and reach an equilibrium 
height, bouncing continues and gradually compactifles the 
aggregates. The enlargement parameter at i = 10"* yr is 
between 2 and 100. The dust distribution reaches a quasi- 
steady state at t = IC* yr. 

We see that the particle mass is constant as a function 
of height at t = 10"* yr. As turbulence effectively mixes the 
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particles, and as bouncing prevents further growth or frag- A value of a = 10 ^ results in height-dependent particle 
mentation (the dust growth is halted), both the masses and mass. 



porosities of the aggregates are similar at all heights where 



We also see from these simulations that a higher turbu- 



dust is present. We will see m the next section that this , i j i.i r i- i j • xi 

, , . r , , , , , , . , , , , ience value reduces the mass of particles and increases the 

IS only true it the turbulence parameter is rather modest. -,, \ \, ■ \,i. re >- -u ^ • >. uiuji 

■^ ^ dust scale height. It turbulence is strong enough, the dust 
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scale height can be similar to the gas scale height and the 
disk atmosphere remains dusty at all times. However, such 
high turbulence value prevents any significant dust growth, 
which is not a fertile environment for planet formation (see 
next Section). 



4.2. The complete Braunschweig collision model 

In this Section we use the complete Braunschweig model 
(see Paper I for details) , the value of the turbulence param- 
eter is a = 10~^ and 10~^. The calculations are performed 
with both the Okuzumi porosity model (CBl, CBS) and 
the Ormel porosity model (CB2, CB4) because we want to 
investigate how sensitive the outcome of dust growth is to 
the used hit&stick porosity model within the context of our 
model. 

In the simplified Braunschweig collision model, the 
growth is halted by bouncing immediately if the parti- 
cles enter the bouncing regime. However, in the complete 
Braunschweig model, there channels for growth beyond the 
hit&stick border line (that is where Ecow > S-Eioii)- The 



most important area is in the 'pP' regime where a small 
porous projectile collides with a heavy porous target (see 
Fig. 11 of Paper I). Due to these "green" areas at inter- 
mediate collision energies, particles in the CBl and CB2 
simulations grow to higher masses than in the SB3 simula- 
tion. As a consequence, the scale height of the dust is lower 
in these simulations, as heavier particles are more difficult 
to stir up by turbulence. We illustrate the mass distribu- 
tion at t = 1, 10, 100, 316, 10^, 10" yr in Fig.[7]for the CBl 
model. 

The particle evolution has two phases in these simu- 
lations. The first 1000 yr are identical for the SB3 and 
CB1/CB2 simulations, respectively (see also the first five 
snapshots of Figs. [5] and [7]). In this phase, particles start 
sedimenting, and the rain-out particles reach the midplane. 
The dust evolution in the SB3 model halts at this point 
as only bouncing collisions happen. However, during the 
second phase of the CBl and CB2 simulations, particles 
can grow to higher masses because there are areas in the 
parameter space that is favorable for growth somewhat be- 
yond the bouncing barrier (see Paper I and Paper II for a 
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(CBS). The axes and contours are the same as in Fig. [5l 



detailed explanation). Due to this growth, particles reach 
10~^ g in mass for the CBl and CB2 simulations. 

The time evolution of the enlargement parameter is 
quite different in the two cases. Fluffy aggregates are pro- 
duced with enlargement parameter ^ — 10^ when the 
Okuzumi porosity model is used. However, these aggregates 
are strongly compacted by bouncing and their final enlarge- 
ment parameter is between 2 and 20. When the Ormel 
porosity model is considered, the enlargement parameter 
never exceeds 20 and by the end of the simulation the en- 
largement parameter is between 2 and 10. As we see, the en- 
largement parameter evolved through very different ways, 
but the final enlargement parameters are not so much dif- 
ferent for the CBl and CB2 simulations. 

Figure [8] illustrates the collision history of the CBl sim- 
ulation. If we compare this figure to Figs. 5, 7 and 10 of 
Paper II, we see that the features are more smeared out in 
Fig.|8]than in the other figures. As we simulate here several 
boxes at different heights above the midplane, the physical 
conditions (e.g. gas density and sedimentation velocity) at 
the midplane and at the upper scale heights of the disk are 
different, which is responsible for the smeared out features 
of Fig. E 



We investigate the collision frequency of the nine colli- 
sion types as a function of time (see Fig.|9]). If one compares 
this figure with Figs. 4c, 6c, and 9c of Paper II, we imme- 
diately see that the diversity of occurring collision types is 
much greater in the CBl simulation, although the strength 
of the turbulence is the same in all cases (a = lO""*). This 
can be explained by the presence of sedimentation. The par- 
ticle population in a given box is not fixed as in Paper II. 
During the rain-out process, "heavy" particles coming from 
z = 1-2 Hg 'hit' the particle population at the midplane. 
Due to this process, the relative velocity between the small, 
midplane particle population and the generally larger rain- 
out population is increased. Thus collision types can occur 
that require larger collision energies. 

We explore the effects of strong turbulence (CBS and 
CB4 runs). In general we find that a turbulence parame- 
ter of a = 10~^ is able to keep the upper layers of the 
disk atmosphere dusty (see Fig. ITU)) . However the size of 
the particles at the midplane is several order of magnitude 
smaller than in the CBl or CB2 simulations. This is not fa- 
vorable for planetesimal formation via s elf-gravity assisted 
planetesimal forma tion (jjohansen et al.L 2007; Cuzz i et all 
120081: lYoudinl . 120 111) . The high level of turbulence lowers the 
enlargement parameter in the CBS simulation (see Fig. fTTj) . 
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The value of the enlargement parameter does not exceed 
10'^ during the simulation. The final enlargement parame- 
ter is between 2 and 10 for both the CBS and 4 simulations. 

Interestingly, the particle masses as function of height 
are constant in the simulations with a = 10~^ or 10~^, but 
for higher a values, this is not the case. The particles in the 
upper layers are significantly smaller than particles at the 
midplane. The particles are mass-sorted (i.e. heavy parti- 
cles are mostly located at the midplane and small particles 
at the upper layers). As 100 micron-sized particles from 
the midplane is mixed upwards by the strong turbulence, 
it looses mass along the way due to fragmenting collisions. 
Similarly, as a particle moves closer to the midplane, it 
is growing in mass due to sticking collisions. This can be 
explained by the decreasing gas density as a function of 
height. At low densities the stopping time (and therefore 
the relative velocity) of a particle is high. At high den- 
sities the particles are stronger coupled to the gas, their 
stopping time decreases. Therefore the relative velocity be- 
tween these particles are low and they can coagulate. Figure 
[1^ illustrates this effect as fragmentation (Fl, F2 and FS) 
happens frequently at z=2 Hg above the midplane (green 



shade for Fl collisions) but it is rare at the midplane (dark 
blue shade for Fl collisions). 

We can also verify this behavior by comparing the vis- 
cous and collision timescales in this model. The viscous 
timescale is 



iy 






using the parameters of the disk model we get that ty 
yr. The collision timescale can be calculated as 



tr 



1 



Avan' 



(22) 



22 



(23) 



where Av is the relative velocity of the two particles (in this 
case it is 1 m/s, the critical velocity for fragmentation), a 
is the cross section of the two particles (in this case two 
monomers with 1 micron size as such particles are present 
at 4 Hg), n is the number density of the particles (calculated 
at the height of 4 Hg assuming monomer-sized particles). 
Using these values, we get that ic = 8 yr. Therefore it is 
further verified that particles are fragmented while they are 
mixed upwards by turbulence as the timescale for collisions 
is shorter than the timescale for mixing. 
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It is also important to note that the particles at the 
upper layers of the disk have masses of 10"^" g and below 
(size of ~ 1-10 micron). Although the upper layers of the 
disk are not well resolved, the obtained particle sizes are 
comparable to th e sizes required to expl ain the 10 micron 
feature of SEDs (|van Boekel et all . |2003() and this vertical 
dust distribution has reached steady state after t = 2000 
yr. Therefore, we conclude that in the framework of our 
models, high values of turbulence (therefore energetic colli- 
sions with strong bouncing and fragmentation) are needed 
to explain why the disk atmospheres are dusty throughout 
the lifetime of the disk. 



4.3. Different disk models 

We performed simulations in two additional disk models 
to explore the effects of the gas surface density. These are 
the low density (with Eg — 10 g/cm^) and high density 
simulations (with E^ = 1000 g/cm^). See also Tab.[l] sim- 
ulations LDl-2 and HDl-2. The dust to gas ratio, the gas 
scale height and the gas temperature are the same as in all 



previous simulations, we change the turbulence parameter 
only. 

We find that the final mass and Stokes number of the 
particles depend on the gas density (see Tab.[T]). The higher 
gas density increases the particle masses but decreases the 
Stokes numbers. Thus the low density model produces the 
smallest particles (lO"** g if a = lO""*) but these particles 
have the highest Stokes number amongst the three gas den- 
sities (almost 10~^ if a = lO"**). This can be explained by 
considering the stopping time of these particles. The stop- 
ping time is inversely proportional to the gas density (see 
Eq. [TU| . Therefore particles in the low density model have 
greater stopping times than particles in the high density 
model. Particles with large stopping times are loosely cou- 
pled to the gas and have therefore high collision velocities. 
For this reason, the maximum particle mass proportional 
to the gas density. 

The common feature of these simulations is that the 
dust scale height for a ~ 10^^ is always similar to the gas 
scale height (see Table [ij. Therefore, we conclude that the 
disk atmospheres can be kept dusty in sparse and dense 
disks alike with sufficiently high turbulence values. 
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5. Discussion and future work 

5.1. Validity of tlie porosity model 

The Braunschweig coUision model is based on laboratory 
experiments that used fractal dimension 3 aggregates with 
an enlargement parameter typically between 7 and 2. As 
seen in Table [TJ the maximum enlargement parameter of the 
aggregates can be several orders of magnitude larger than 
the dust aggregates used in the laboratory. More specifi- 
cally, the hit&stick collision regime is followed by a collision 
type called sticking through surface effects (see Paper I) 
and bouncing at higher collisional energies. Thus we assume 
that the transition from fractal aggregates in the hit&stick 
regime to fractal dimension 3 aggregates in bouncing is 
an instantaneous one. Furthermore, our porosity model for 
bouncing is calibrated for dust cakes (enlargement parame- 
ter of 6), however the enlargement parameter at the end of 
the hit&stick phase is 2-3 orders of magnitude larger than 
the enlargement parameter of the dust cakes (see Table [T]). 
Therefore, it is questionable whether the porosity model for 
the restructuring phase is correct. Significant changes in the 
porosity of aggregates has the potential to significantly al- 
ter our collision model and thus the obtained results. 

It is still debated how the porosity and the fractal di- 
mension evolves during the evolution of dust aggregates. As 
this is a central issue in determining the dust properties, we 
need detailed information on the porosity evolution of dust 
aggregates. To resolve these issues, we plan to follow the 
hit&stick and restructuring phases of a particle distribu- 
tion in a self-consistent way by combining the Monte Carlo 
model of ZsD08 with a molecular dynamics code in a follow- 
up paper. 



5.2. Aggregates beyond the snow line 

The particle sizes at the midplane of the disk are rather 
small (~ 100 microns if a = 10~^) as we see in the pre- 
vious section. The question arises: under what conditions 
could planetesimals form via successive collisions of dust 
aggregates? Or how do large enough dust aggregates form 
that could, under favorable conditions, be concentrated in 
e.g. vortices or turbulent eddies, become self-gravitating, 
and form ev entua lly planetesimals (Johansen et al., 2 007; 
ICuzzi et all . 120081: lYoudinl . 120111: iJohansen et all . 1201 Ih ? 

One answer to these questions coul d be icy aggregates . 
The molecular dynam ics simulations of lWada et al.l (|2008D ; 
ISuvama et al.l (|2008[) showed that icy aggregates are very 
resilient to restructuring. They observed sticking between 
icy aggregates at a relative velocity as high as 50 m/s. 
The uncertainty in these simulations is the microphysi- 
cal parameters of the icy monomers (such as critical dis- 
placement, surface energy. Young modulus etc.). Recently 
laboratory experiments were performed using micron-sized 
ice monomers by [Gundlach et al. (2011J). They measured 
the rolling energy between icy mono mers and it turne d 
out the previously a ssumed values by IWada et all (|2008t) ; 
ISuvama et al.l ()2008[ ) are in good agreement with the mea- 
sured laboratory values. If experiments also confirm that 
icy aggregates can stick at relative velocities as high as 50 
m/s, that would provide a way to form large enough dust 
aggregates beyond the snow-line in the solar nebula. 



5.3. Is a constant as a function of height? 

iGammiel ()1996[) proposed the concept of layered accre- 
tion disks. If the ionization fraction of the gas is not suf- 
ficient to support magneto-rotational instability (MRI - 
iBalbus & HawlevI (|l99ll) ). the turbulence parameter drops 
and a dead-zone forms at the midplane of the disk. The 
extent of the dead-zone is uncertain, as the ionization 
processes of the gas are not well-constrained. For typ- 
i cal TTauri disks it can extend between 0.1-4 AU 
(jD'Alessio et all . I1998D . Inside 0.1 AU, the thermal radi- 
ation from the star can keep the gas sufficiently ionized for 
MRI, and outside 4 AU, the gas surface density is typically 
below 100 g/cm^, therefore cosmic rays can penetrate the 
disk and keep it MRI acti ve at all h e ights. 

It was proposed by lOkuzumil (J2009t) that negative 
charges on the surfaces of grains could prohibit growth. As 
we use a Monte Carlo code to follow the evolution of aggre- 
gates, it is possible to include a third particle property: the 
amount of charges present on the grain. In order to follow 
the charge evolution of the grains, we need to solve for the 
ionization state of the gas (i.e. amount of charges available 
in the gas phase), how efficiently these charges are collected 
by the dust grains and how the charges affect the relative 
velocity of the aggregates. 

We plan to investigate how dust evolves in a layered 
disk model. Small dust particles can very efficientl y swee p 
up charges in the gas. As shown bv I Turner et al.l (|2010t ). 
the dead-zone can extend to 2 Hg for 1 micron-sized parti- 
cles, but it shrinks below 0.5 Hg for aggregates that are 100 
micron in size. In a simulation like the one presented here, 
this would mean that as the particles grow, the dead-zone 
shrinks. When the dead-zone disappears, the whole disk be- 
comes MRI active and the particles settled to the midplane 
might be fragmented and stirred back up. This could lead 
to initial oscillations before an equilibrium state is reached. 

6. Summary 

We performed simulations in a ID vertical column of a pro- 
toplanetary disk to better understand the process of sedi- 
mentation. We simultaneously solved for the particle mo- 
tion and growth inside this column. The complexity of the 
models was gradually increased to examine the effects of dif- 
ferent processes. The first simulation used a collision model 
that only contained sticking. We furthermore assumed that 
the particles were compact, and the turbulence parame- 
ter (a) was set to zero. Later on we investigated the ef- 
fects of different porosity models, more realistic collision 
models (with sticking, bouncing and fragmentation) and 
turbulence of different strengths. Below we summarize our 
results. 

— Porosity helps to produce heavier particles by decreas- 
ing the stopping time of the aggregates, but porosity 
alone cannot prevent rain-out. 

— If the size of the particle becomes greater than the mean 
free path of the gas, the drag law changes from Epstein 
to Stokes drag. At this point, rain-out becomes unavoid- 
able due to the change in the drag law. 

— Realistic collision models with bouncing and fragmen- 
tation limit the maximum particle sizes to be not more 
than 1 mm - 1 cm (particle masses between 10~^ - 1 g 
). The exact value depends strongly on the strength of 
the turbulence, and the gas density. 
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— A higher value of turbulence decreases the particle 
masses but increases the dust scale height. Using the 
simplified Braunschweig model with a — fO^® results 
in a dust scale height of 0.2 Hg and final particle mass 
of 10~^ g, however the dust scale height is 0.6Hg, and 
the final particle mass is fO~^ g when using a = 10"''. 

— The final particle size and Stokes number depends on 
the gas density. The mass of the particles is decreasing 
with decreasing gas density, however the Stokes number 
remains roughly constant. 

— When using the most detailed collision model (the com- 
plete Braunschweig collision model), we obtain particle 
masses of 10~^ g (with an average radius of 1 mm, and 
an average Stokes number of 4 x 10""^) and a dust scale 
height of 0.2 Hg upon using a = 10""*. However, the 
dust scale-height is almost 1 Hg and the final particle 
mass at the midplane is 10~^ g (with an average radius 
of 100 micron, and an average Stokes number of 10~^) 
upon using a = 10"^. Therefore, a sufficiently high tur- 
bulence value can keep the disk atmosphere dusty but 
the absence of significant dust growth is not favorable 
for planet formation. 

— The maximum enlargement parameter of the aggregates 
during their evolution can be as high as 10*. However, 
our adopted porosity model for bouncing is based on 
laboratory experiments performed using enlargement 
parameter ~ 6 dust cakes. Therefore the largest uncer- 
tainty of our model is the porosity evolution. 

— The dust particle mass as a function of height is not 
constant if the turbulence parameter is a = 10~^. We 
see that the particle mass/size is a decreasing function 
of the height. Particles are 100 microns in size at the 
midplane and a few microns at 4 pressure scale heights. 

— The micron-sized particles present in the upper layers 
are comparable to the sizes needed to explain the 10 
micron feature of disk SEDs. Therefore in the framework 
of our model, high values of turbulence are needed to 
explain why disk atmospheres are dusty for ^ 10^ yr. 
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